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Abstract -We present an analytical closed form expression, which gives a good approximate 
propagator for diffusion on the sphere. Our formula is the spherical counterpart of the Gaussian 
propagator for diffusion on the plane. While the analytical formula is derived using saddle point 
methods for short times, it works well even for intermediate times. Our formula goes beyond 
conventional "short time heat kernel expansions" in that it is nonperturbative in the spatial 
coordinate, a feature that is ideal for studying large deviations. Our work suggests a new and 
efficient algorithm for numerical integration of the diffusion equation on a sphere. We perform 
Monte Carlo simulations to compare the numerical efficiency of the new algorithm with the older 
Gaussian one. 



' O Introduction: Diffusion on the sphere is a problem 
^ that arises in several contexts. At the cellular level, dif- 

o 

^fusion is an important mode of transport of substances, 
i— The cell wall is a lipid membrane and biological substances 
like lipids and proteins diffuse on it. In general biological 
^.membranes are curved surfaces. Spherical diffusion also 
QO crops up in the swimming of bacteria, surface smoothen- 
f* — ing in computer graphics [l] and global migration patterns 
O^l of marine mammals [2] . Such diffusion effects are stud- 
.ied using computer simulations. There have been experi- 
C*~) mental studies of fluorescent marker molecules on curved 
l O surfaces like micelles and vesicles [3|[4j using fluorescence 
anisotropy decay. 

J* While diffusion in Euclidean space has been studied ex- 
. intensively both analytically and numerically [5}|6], there 
^have been fewer studies of diffusion on curved surfaces. 
There have been some analytical studies of diffusion on 

03 



[7J|8] arrive at 
for the proba- 



curved surfaces [7^(9] . Some of these studies 
a short time heat kernel expansion 10 
bility distribution. This expansion is perturbative in the 
time step as well as the spatial step as is evident from 
Eq. (51) of Ref. [8]. However, there does not exist a sim- 
ple closed form analytical expression, for the probability 
distribution for diffusion on a curved surface which can 
be implemented in a computer simulation. Here we fo- 
cus on diffusion on a sphere, the case of a curved surface 



with constant positive curvature. The simplest and most 
natural method of generating diffusion on a curved man- 
ifold is to consider a large number of walkers performing 
a random motion. The random walkers can be described 
by a Langevin equation using a spatial step size much 
smaller than the inverse curvature so that curvature can 
be neglected. In this "tangent space" approximation, the 
sphere appears locally flat and one can ignore its curva- 
ture. This is the approach taken for example in [3j|4] in 
which the diffusing particle takes a Brownian step in the 
tangent plane at a point on S 2 and is then projected back 
to the surface of the sphere. One can describe the ensem- 
ble of random walkers with a probability distribution P 
and write a Fokker-Planck (diffusion) equation governing 
the spread of probability. One finds (see page 847 of Ref. 



13 ) on general grounds, the temporal step size must be 



chosen smaller than the square of the spatial step size. As 
a result of the small spatial step size forced on us by the 
tangent space method, integrating in time is computation- 
ally expensive. It would be useful to have a method which 
permits larger step sizes, not constrained by the curvature. 



There has been some work 14.15 generating an algo- 



rithm for diffusion on S 2 by making use of known ana- 
lytic results on S 3 , the hypersphere in four dimensions 
and identifying S 2 with a thin strip around the equator 
of S 3 with a reflecting boundary. These methods give 
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rise to an improved algorithm for treating diffusion on S 2 . 
However, the use of S 3 and the method of dimensional 
reduction render the treatment less than intuitive. Our 
purpose here is to provide a more down to earth approach 
by working directly on S 2 . We derive an approximate an- 
alytical formula for the distribution function representing 
diffusion on a sphere. In contrast to earlier studies [7j[8] 
we use a semiclassical saddle point approximation which 
is perturbative in time but nonperturbative in space. Our 
formula reduces to the formula presented in [7j[8j in the 
limit of small spatial scales. In the planar limit our for- 
mula reduces to a Gaussian. We use our "Gaussian" on 
the sphere to generate an efficient algorithm for simulat- 
ing Brownian motion on a sphere. The algorithm improves 
over the earlier tangent space simulation algorithm. We 
expect this algorithm to have a wide range of applications, 
especially in fields where one needs to repeatedly generate 
diffusion processes on spheres. In the physics of stiff, inex- 
tensible polymers, one has contraints on the bond lengths 
while the bond angles are relatively unconstrained. This 
leads to random motions, as in the Kratky-Porod model 
which represent diffusion on spheres. In computer graph- 
ics, one seeks algorithms which can be used to smooth 
graphical data on spherical images. A standard technique 
is to use the diffusion equation on the sphere, but a simple 
and accurate algorithm for implementing this is presently 
lacking. 

Theory: To define our problem regarding diffusion on 
a sphere, consider first diffusion on a plane which is de- 
scribed by 



dP 
~dt 



(x,t) = DV 2 P{x,t) 



(1) 



where D is the diffusion constant and P(x, t) the proba- 
bility distribution on the plane at time t. The elementary 
solution with initial distribution P(x, 0) = S 2 (x) is easily 
found by Fourier transform. The formal solution is 



P(x,t) 



1 

4^ 



dk e~ Dk - kt e ik - ] 



(2) 



which can be integrated to give a simple Gaussian form 

(3) 



, . / 1 \ x • x 

P&t) = ( 737 ^ ] exp- 



ADt 



When expressed in plane polar co-ordinates (r, 9), the 
probability distribution function in r is 



/ T \ T 

P( ^ ) = (2^) eXP 



4Dt 



(4) 



While these two expressions (2j4| are formally equal to 
each other, the first is an unwieldy integral whereas the 
second is a simple closed form Gaussian expression which 
is well suited to numerical work. In ([2| the time t appears 
in the numerator of the exponent leading to slow conver- 
gence for small times, unlike Q where it appears in the 
denominator. Our purpose in this Letter is to seek the 
analogue of Q for spherical diffusion. 



Consider the analogous situation on the sphere 

S 2 = {(x\x 2 ,x 3 ) eM 3 : (x 1 ) 2 + (x 2 ) 2 + (x 3 ) 2 = R 2 } 

(5) 

of radius R in three dimensional Euclidean space. In stan- 
dard polar coordinates on the sphere the diffusion equation 

dP_D_ (J_^_ ■ n<^ 1 d 2 P \ 

dt ~ r 2 {^8W sm ~m + ^Te W) 

Let us define a dimensionless time variable r = 2Dt/R 2 
for convenience. Typical values for diffusion on a lipid 
membrane are R — 1/im and D = l/xm 2 /sec. The r values 
of interest to us are on the order of seconds. 

For an initial 9 distribution P(9, 0) (including the mea- 
sure sin 9) with S function support at the North pole 9 = 0, 
the diffusion equation has a formal solution for the final 
distribution in 9 at time r: 



P(6,t) 



sin 9 v ^ 

1=0 



{2l + l)e- l(l+1 ^ T / 2 P l (cos9). (7) 



Such a solution is formally exact and analogous to 
but in practice, unwieldy to use in numerical work since 
it is expressed as an infinite series. For short times, the 
convergence of the series is poor and truncation leads to 
spurious oscillations. 

Is there an analogue of the Gaussian form Q on the 
sphere? The main result of this Letter is a closed form ap- 
proximate expression which generalises the planar Gaus- 
sian Q to the sphere. This approximate expression is 
given by: 

Q{9,T) = Af ^V(h^8e-£ (8) 

T 

where A/"(r) is a constant determined by the normalisation 
condition 

•) = 1 (9) 



In the limit of small times, the expression ([8]) reduces to 
i#exp — ^ which agrees with the planar Gaussian 
Our study shows that this expression turns out to be a 
good approximation to the exact propagator ([7]) even for 
intermediate times. 

The Heat Kernel Expansion (HKE) results in a propa- 
gator (see Eq. 51 in Ref.[ M]) expressed as a power series 
in 9 2 

Qhke(0,t) = -(1 + — + .. .)e~£ sin# (10) 

T LZ 

where we have transcribed the content of Eq. 51 in Ref.[ 
[8]] (suitably restricted to constant curvature) and inserted 
the measure factor of sin 9 for easy comparison. The new 
propagator is essentially a summation of this perturbative 
series. 

Figures 1 and 2 show a comparison of the exact prop- 
agator (with the upper limit l m ax in the summation set 
to 20) with the approximate one for a range of t. The 
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plots of Fig. show P(9, t) and Q(9, r) for r = .5, 1 
and Q{9,t)hke truncated to the second term. 

The deviation between probability distributions can be 
quantitatively measured using the Kullback-Leibler diver- 
gence, which is also known as the relative entropy. This 
measure which is widely used in information theory |16| 
gives a positive number Dkl which measures the extent 
of deviation (or divergence) between a trial distribution 
Q(9) and a fiducial one P{9): 



Dkl := 



(11) 



Dkl vanishes if and only if the two normalised distribu- 
tions are identical. Dkl is also invariant under changes 
of coordinates. For the fiducial distribution we use the ex- 
act form ([7]) and compare KL divergence of the Gaussian 
propagator (restricted to the interval [0, tt] and normalised 
in this interval), the Qhke (10) as well as our new dis- 
tribution Qnew The table shows a comparison of Dkl 
Gaussian, Dkl HKE and Dkl new, for different values 

Of T. 

Note that the new propagator always has a smaller KL 
divergence than the Gaussian propagator and Qhke ■ It 
is thus a better approximation to the exact propagator. 



time t 


Dkl Gaussian 


Dkl HKE 


D K l new 


0.1 
0.5 
0.7 
0.9 
1.0 
2.0 


1.4x 10~ 4 
4.0x 10~ 3 
8.2x 10" 3 
1.3x 10" 2 
1.6x 10~ 2 
4.3x 10~ 2 


6.5 x 10~ 7 
3.7 x 10- 4 
1.3x 10~ 3 
2.9x 10~ 3 
3.9 x 10" 3 
l.lx 10~ 2 


1.8 x 10~ y 
3.8 x 10" 6 
2.9x 10~ 5 
9.9x 10~ 5 
1.6 x 10~ 4 
8.9x 10~ 4 



Derivation of the Main Result: The elementary so- 
lution to the diffusion equation on the unit sphere is given 



by the Wiener integral 17 



K{n u 0, n 2 , r) = fv[h(T)]exp-S[n(T)} (12) 



where 



S[h(r)} 



dh dh 
dr dr 



-dr. 



(13) 



K(hx, 0, h 2l t) is the conditional probability that the 
particle will be at h 2 at time r given that it was at hi 
at time 0. For short times r , we may use a semiclassical 
approximation |18| 



K(hi, 0, h 2 , t) = Af(r) VdetV exp-S, 



ci \ni,n 2 ,T\ 



(14) 



where S c i[hi,h 2 ,T)] is defined as the classical action of 
the least action path connecting hi to hi in "time" r. 
The least action path is unique if hi,h 2 are not collinear, 
which we assume and thus exclude the isolated points 9 = 
0, tt. Det V is the Van Vleck determinant given by the 
determinant of the 2x2 Hessian matrix 



d 2 S, 



i n 2 , 



dh\ dh 2 



(15) 



and Af is a normalisation constant. 

Varying the action to find the classical path yields the 
geodesic equation 

d 2 h 



= An, 



(16) 



where A is a Lagrange multiplier enforcing the constraint 
h ■ h — 1. The solution to (16 1 is the unique great circle 



passing through hi and h 2 . The classical action is given 
by 

5* cl = 9 2 /2t (17) 

where cos 9 — hi ■ h 2 defines 9, < 9 < tt, the length of 
the shortest geodesic arc connecting hi to h 2 . From an 
evaluation of the Van Vleck determinant 



detV 



hip e pd h 2q 



- qok Vu V lk 



wc find that 



detV 



t 2 sin 9 

This leads to the approximate propagator 

Af(r) 



K[hi,Q,h 2 ,r] 



smf e 2t 



(18) 



(19) 



Multiplying by the measure sin 9 to convert into a 9 distri- 
bution leads to the approximate propagator which is used 
in eq.([8|. 

Computer Simulations: We have performed Monte 
Carlo simulations of spherical diffusion to investigate the 
numerical efficacy of the approximate propagator derived 
analytically in the last section. The Fokker-Planck equa- 
tion on a sphere is exactly solvable and the solution can be 
expressed as a series involving the Legendre polynomials as 
described earlier. However, because of poor convergence 
of the series for short times, this form is not suitable for 
repeated numerical evaluation. We use the exact solution 
with a comfortably large cutoff l ma x — 20 as a standard 
against which the computer simulations are tested. 

In performing the simulations, we have closely followed 
the method adopted in Ref. [3j|4] for studying diffusion 
of fluorescent molecules on a sphere. We have studied the 
diffusion of 5 x 10 6 molecules and obtained the distribution 
of the final polar angular displacement 9. 

To integrate for a time r, we split up the time interval 
r into n step time steps t ~ n stev T s t ep each of size r step . 
The value of T s t ep determines the spatial step size a via 
the standard diffusion relation 



'step 



a 

Y 



(20) 



We have implemented two distinct algorithms for ran- 
domly choosing f3: a Gaussian algorithm and our pro- 
posed new algorithm. In the first algorithm, we replace 
the spherical geometry locally by the tangent plane to 
the sphere at the starting point of each Monte Carlo 
move. Confining ourselves to the local tangent plane on 
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the sphere at each move, we choose j3 according to the 
distribution (eq. |3| 



a variable coordination number on the lattice to encode 
curvature 



20 



Q 



Gauss 



(0) 



2/3 



exp(-^> 2 )- 



(21) 



Thus at each step, the angular displacement /3 is chosen 
from a two dimensional Gaussian distribution which has 
zero mean and standard deviation a. This is our first 
algorithm. It is only expected to be accurate when the 
step size a is small. 

In the second method, the angular displacement at each 
step was chosen from the approximate propagator given by 

©• 

2J\f 

Qnem{P) = V SHI f3 exp (- /3 2 / CT 2 ) (22) 



This distribution reduces to the previous one (21) when 



the step size a is small, but has a larger range of validity 
since it takes into account the curvature of the sphere. 

In both algorithms the random numbers were generated 
by constructing a random walk of the variable (3 in the ex- 
ternal potential —logQ(/3). After discarding 10 5 steps to 
get rid of transients, we arrange the subsequent 10 6 j3 val- 
ues as a 10 3 x 10 3 matrix, transpose it and store the values 
for use in the program. This shuffles the random numbers 
and removes correlations between neighboring steps. 

Whereas the planar Gaussian algorithm needs to use 
small displacements to ensure that the tangent plane ap- 
proximation remains valid, the proposed new algorithm is 
not constrained by this requirement. It remains a good 
approximation even for intermediate a. As a result in or- 
der to integrate for r = 2, we need to use just two time 
steps with a step size of a = 1.414 with the new algorithm 
(Figure 3). With the Gaussian algorithm, even with eight 
time steps of step size a — .7071, one achieves a poorer 
accuracy (Figure 4). This new algorithm is our second 
main point in this Letter. 

Conclusion. — In this Letter we derive an approxi- 
mate analytical formula for the distribution function for 
spherical diffusion and use it to generate a simple and 
efficient algorithm for simulations. While we have re- 
stricted ourselves to spherical diffusion, it is evident that 
our results apply equally well to the saddle, (or hyperbolic 
plane), the space of constant negative curvature. The only 
change is that the circular function sin0 is replaced by 
the hyperbolic function sinh#. Our results for the ap- 
proximate propagator ( 19 ) also extend to higher (D) di- 



mensional spheres and saddles with the slight modification 
that the prefactor #/sin# (or its hyperbolic form) is re- 
placed by (8/ sintf)'^. 

Rasin et al [19] studied diffusion on a planar lattice 
and addressed the problem of scaling of the time step as 
the square of the lattice size. Their solution uses kinetic 
methods that replace the parabolic diffusion equation by 
a hyperbolic one, thus allowing for larger time evolution 
steps. It would be interesting to extend their work, using 



For a general curved surface embedded in three dimen- 
sional space, the intrinsic geometry is determined entirely 
by the Gaussian curvature k. Since diffusion on the sur- 
face depends only on the intrinsic geometry, we may con- 
sider each part of the surface as locally approximated by 
a sphere of radius R — We may then apply the 

new simulation algorithm described in this paper in a lo- 
cal patch. This entails numerically introducing geodesic 
based coordinates and choosing isotropically distributed 
steps, with a size distribution given by (22). We expect 



this to improve on the tangent space simulation methods, 
since these only approximate the local geometry to first 
order and do not take into account the second order ef- 
fects of curvature. The new propagator would be useful in 
a wide variety of applications. 



It is a pleasure to acknowledge many stimulating dis- 
cussions with Abhishek Dhar, Rajaram Nityananda and 
Sanjib Sabhapandit. 
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Fig. 1: (Colour online) Figure shows a comparison of the exact propagator (Eq. |7| truncated to l max = 20), the approximate 
one Eq. |8} and the Heat Kernel expansion for r = 0.5, 1.0. Note that the difference between the first two is barely discernible. 




Fig. 2: (Colour online) Figure shows a comparison of the exact propagator (Eq. |7| truncated to l ma x = 20) with the approximate 
one (Eq. Q) for r = 1.5. The difference is now apparent. Also shown is the Heat Kernel expansion, which differs even more 
from the exact propagator. 



p-6 



A "Gaussian" for Diffusion on the sphere 




Fig. 3: (Colour online) Figure shows a comparison of the exact propagator with the simulation data using the new algorithm. 
The total time r of integration is r = 2 and this time has been achieved with two steps, each of step size a = 1.414. 




Fig. 4: (Colour online) Figure shows a comparison of the exact propagator with the simulation data using the Gaussian 
algorithm. The total time r of integration is r = 2 and this time has been achieved with two steps of step size a = 1.414 
(dashed curve) and eight steps of size a = .7071. Notice that the Gaussian algorithm requires a larger number of steps to 
reproduce the exact analytical propagator. 
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